home *** CD-ROM | disk | FTP | other *** search
/ Technotools / Technotools (Chestnut CD-ROM)(1993).ISO / lang_c / cug167 / spline.c < prev    next >
Text File  |  1985-08-21  |  15KB  |  658 lines
  1. /* 
  2. HEADER:     CUG
  3. TITLE:        SPLINE.C - Interpolate Smooth Curve
  4. VERSION:    1.00
  5. DATE:        11/01/85
  6. DESCRIPTION:    "SPLINE takes pairs of numbers from the standard
  7.         input as abscissae and ordinates of a function.
  8.         It produces a similar set, which is approximately
  9.         equally spaced and includes the input set, on the
  10.         standard output. The cubic spline output has two
  11.         continuous derivatives and sufficiently many
  12.         points to look smooth when plotted."
  13. KEYWORDS:    spline, graphics, plot, filter, UNIX
  14. SYSTEM:        ANY
  15. FILENAME:    SPLINE.DOC
  16. WARNINGS:    NONE
  17. CRC:        xxxx
  18. SEE-ALSO:    SPLINE.C
  19. AUTHORS:    Ian Ashdown - byHeart Software
  20. COMPILERS:    ANY
  21. REFERENCES:    AUTHORS: Bell Laboratories
  22.         TITLE:     UNIX Programmer's Manual, Volume 1
  23.              p. 145; Holt, Rinehart and Winston
  24.         AUTHORS: R.W. Hamming
  25.         TITLE:     Numerical Methods for Scientists and
  26.              Engineers, 2nd Edition
  27.              pp. 349 - 355; McGraw-Hill (1973)
  28.         AUTHORS: Forsythe, Malcolm & Moler
  29.         TITLE:     Computer Methods for Mathematical
  30.              Computations
  31.              pp. 70 - 83; Prentice-Hall
  32. ENDREF
  33. */
  34.  
  35. /*-------------------------------------------------------------*/
  36.  
  37. /* SPLINE.C - Interpolate Smooth Curve
  38.  *
  39.  * Version 1.00        November 1st, 1985
  40.  *
  41.  * Copyright 1985:    Ian Ashdown
  42.  *            byHeart Software
  43.  *            1089 West 21st Street
  44.  *            North Vancouver, B.C. V7P 2C6
  45.  *            Canada
  46.  *
  47.  * This program may be copied for personal, non-commercial use
  48.  * only, provided that the above copyright notice is included in
  49.  * all copies of the source code. Copying for any other use
  50.  * without previously obtaining the written permission of the
  51.  * author is prohibited.
  52.  *
  53.  * Synopsis:    SPLINE [option] ...
  54.  *
  55.  * Description:    SPLINE takes pairs of numbers from the standard
  56.  *        input as abscissae and ordinates of a function.
  57.  *        It produces a similar set, which is approximately
  58.  *        equally spaced and includes the input set, on the
  59.  *        standard output. The cubic spline output (R.W.
  60.  *        Hamming, "Numerical Methods for Scientists and
  61.  *        Engineers", 2nd ed. 349ff) has two continuous
  62.  *        derivatives and sufficiently many points to look
  63.  *        smooth when plotted.
  64.  *
  65.  *        The following options are recognized, each as a
  66.  *        separate argument:
  67.  *
  68.  *        -a  Supply abscissae automatically (they are
  69.  *            missing from the input); spacing is given by
  70.  *            the next argument, or is assumed to be 1 if
  71.  *            next argument is not a number.
  72.  *
  73.  *        -k  The constant "k" is used in the boundary
  74.  *            value computation
  75.  *
  76.  *            y " = ky " , y " = ky "
  77.  *             0    1     n         n-1
  78.  *
  79.  *            is set by the next argument. By default,
  80.  *            k = 0.
  81.  *
  82.  *        -n  Space output points so that approximately n
  83.  *            intervals occur between the lower and upper
  84.  *            "x" limits (default n = 100.)
  85.  *
  86.  *        -p  Make output periodic, i.e. match derivatives
  87.  *            at ends. First and last input values must
  88.  *            agree.
  89.  *
  90.  *        -x  Next 1 (or 2) arguments are lower (and upper)
  91.  *            "x" limits. Normally these limits are
  92.  *            calculated from the data. Automatic abscissae
  93.  *            start at lower limit (default 0).
  94.  *
  95.  * Diagnostics:    When data is not strictly monotone in "x", SPLINE
  96.  *        reproduces the input without interpolating extra
  97.  *        points.
  98.  *
  99.  * Bugs:    A limit of 1000 input points is silently
  100.  *        enforced. 
  101.  *
  102.  * (Note: the above description is a slightly reworded version of
  103.  *       that appearing in the "UNIX Programmer's Manual",
  104.  *      copyright 1979, 1983 Bell Laboratories.)
  105.  */
  106.  
  107. /*** Definitions ***/
  108.  
  109. #define FALSE         0
  110. #define TRUE         1
  111. #define MAX_SIZE  1000    /* Input point array limit */
  112.  
  113. #define ILL_ARG         0    /* Error codes */
  114. #define ILL_CMB      1
  115. #define ILL_KVL         2
  116. #define ILL_OPT      3
  117. #define INS_INP      4
  118. #define MIS_KVL      5
  119. #define MIS_XVL      6
  120. #define MIS_YVL         7
  121. #define NMT_ORD         8
  122.  
  123. #define SQUARE(a) a*a
  124. #define CUBE(a) a*a*a
  125.  
  126. /*** Typedefs ***/
  127.  
  128. typedef int BOOL;    /* Boolean flag */
  129.  
  130. /*** Include Files ***/
  131.  
  132. #include <stdio.h>
  133. #include <ctype.h>
  134. #include <math.h>
  135.  
  136. /*** Main Body of Program ***/
  137.  
  138. int main(argc,argv)
  139. int argc;
  140. char **argv;
  141. {
  142.   int n = 0,
  143.       i,
  144.       j,
  145.       n_val = 0;
  146.   float    x[MAX_SIZE],
  147.     y[MAX_SIZE];
  148.   double a_val = 1.0,
  149.      k_val = 0.0,
  150.      x1_val = 0.0,
  151.      x2_val = 0.0,
  152.      x_intvl,
  153.      ix,
  154.      iy,
  155.      d2y[MAX_SIZE],
  156.      h,
  157.      atof(),
  158.      fabs(),
  159.      spl_int();
  160.   char buffer[257],
  161.        *temp,
  162.        *gets();
  163.   BOOL aflag = FALSE,    /* Command-line option flags */
  164.        kflag = FALSE,
  165.        pflag = FALSE,
  166.        x1flag = FALSE,
  167.        x2flag = FALSE,
  168.        is_float();
  169.   void spl_coeff(),
  170.        pspl_coeff(),
  171.        error();
  172.  
  173.   /* Parse the command line for user-selected options */
  174.  
  175.   while(--argc)
  176.   {
  177.     temp = *++argv;
  178.     if(*temp != '-')    /* Check for legal option flag */
  179.       error(ILL_OPT,*argv);
  180.     else
  181.       switch(toupper(*++temp))
  182.       {
  183.     case 'A':    /* "-a" option */
  184.       aflag = TRUE;
  185.       if(argc > 1 && is_float(*(argv+1)))
  186.       {
  187.         argc--;
  188.         argv++;
  189.          if((a_val = atof(*argv)) < TINY)
  190.           error(ILL_ARG,*argv);
  191.       }
  192.       break;
  193.     case 'K':    /* "-k" option */
  194.       if(pflag == TRUE)
  195.         error(ILL_CMB,NULL);
  196.       kflag = TRUE;
  197.       if(argc > 1 && is_float(*(argv+1)))
  198.       {
  199.         argc--;
  200.         argv++;
  201.         k_val = atof(*argv);
  202.         if(k_val - 4.0 < 0.000001)
  203.           error(ILL_KVL,*argv);
  204.         break;
  205.       }
  206.       else
  207.         error(MIS_KVL,NULL);
  208.     case 'P':    /* "-p" option */
  209.       if(kflag == TRUE)
  210.         error(ILL_CMB,NULL);
  211.       pflag = TRUE;
  212.       break;
  213.     case 'X':    /* "-x" option */
  214.       x1flag = TRUE;
  215.       if(argc > 1 && is_float(*(argv+1)))
  216.       {
  217.         argc--;
  218.         argv++;
  219.         x1_val = atof(*argv);
  220.       }
  221.       else
  222.         error(MIS_XVL,NULL);
  223.       if(argc > 1 && is_float(*(argv+1)))
  224.       {
  225.         x2flag = TRUE;
  226.         argc--;
  227.         argv++;
  228.         x2_val = atof(*argv);
  229.       }
  230.       break;
  231.     default:    /* "-n" option */
  232.       while(isdigit(*temp))
  233.         n_val = n_val * 10 + (*temp++ - '0');
  234.       if(*temp != '\0')
  235.         error(ILL_OPT,*argv);
  236.       }
  237.   }
  238.   if(n_val == 0)    /* Set "n_val" if not given */
  239.     n_val = 100;
  240.  
  241.   /* Get the input data */
  242.  
  243.   while(1)        /* ... while there is more input data */
  244.   {
  245.     if(aflag == TRUE)    /* Automatic abscissae were called for */
  246.     {
  247.       if(n == 0)
  248.     x[0] = x1_val;
  249.       else
  250.     x[n] = x[n-1] + a_val;
  251.     }
  252.     else        /* Abscissae supplied with input data */
  253.     {
  254.       if(gets(buffer))
  255.     x[n] = atof(buffer);
  256.       else
  257.     break;
  258.     }
  259.     if(gets(buffer))    /* Read in the corresponding ordinate */
  260.       y[n] = atof(buffer);
  261.     else
  262.     {
  263.       if(aflag == TRUE)
  264.     break;
  265.       else
  266.     error(MIS_YVL,NULL);
  267.     }
  268.     if(++n == MAX_SIZE)    /* Maximum amount of input data? */
  269.       break;
  270.   }
  271.   if(n <= 2)        /* Check for insufficient input data */
  272.     error(INS_INP,NULL);
  273.  
  274.   /* Check for non-monotonic abscissae. Output input data set
  275.    * without interpolation if true.
  276.    */
  277.  
  278.   h = x[1] - x[0];    
  279.   for(i = 1; i < n-1; i++)
  280.     if(fabs(x[i+1] - x[i] - h) > 0.0)
  281.     {
  282.       for(i = 0; i < n; i++)
  283.     printf("%g\n%g\n",x[i],y[i]);
  284.       exit();
  285.     }
  286.  
  287.   /* Calculate abscissa interval. Use "-x" option values if
  288.    * they were given. 
  289.    */ 
  290.  
  291.   if(x1flag == FALSE)
  292.     x1_val = x[0];
  293.   if(x2flag == FALSE)
  294.     x2_val = x[n-1];
  295.   x_intvl = (x2_val - x1_val)/(n_val - 1);
  296.  
  297.   /* Find the coefficients */
  298.  
  299.   if(pflag == FALSE)
  300.     spl_coeff(y,d2y,h,n,k_val);
  301.   else
  302.     pspl_coeff(y,d2y,h,n);
  303.  
  304.   /* Interpolate and output results */
  305.  
  306.   ix = x1_val;
  307.   i = 1;
  308.   for(j = 0; j < n_val; j++)
  309.   {
  310.     if(ix >= x[i] && i < n - 1)
  311.       i++;
  312.     iy = spl_int(ix,x,y,d2y,h,i);
  313.     printf("%g\n%g\n",ix,iy);
  314.     ix += x_intvl;
  315.   }
  316. }
  317.  
  318. /*** Functions ***/
  319.  
  320. /* SPL_COEFF() - Calculate spline coefficients and return in
  321.  *         vector "d2y[]". Matrix to be solved has the
  322.  *         typical form:
  323.  *
  324.  *     +-          -+ +-   -+       +-           -+
  325.  *    |  4-k 1   0   0   | | y1" | = m * | y2 - 2*y1 + y0 |
  326.  *    |  1   4   1   0   | | y2" |       | y3 - 2*y2 + y1 |
  327.  *    |  0   1   4   0   | | y3" |       | y4 - 2*y3 + y2 |
  328.  *    |  0   0   1   4-k | | y4" |       | y5 - 2*y4 + y3 |
  329.  *    +-          -+ +-   -+       +-           -+
  330.  *
  331.  *         where k = k_val, m = 6.0/(h*h) and yn" is the
  332.  *         second derivative of the interpolated function
  333.  *         at the "nth" abscissa, y0" = k * y1" and
  334.  *         y5" = k * y4".
  335.  */
  336.  
  337. void spl_coeff(y,d2y,h,n,k_val)
  338. float y[];
  339. double d2y[],
  340.        h,
  341.        k_val;
  342. int n;
  343. {
  344.   double m,
  345.      a[MAX_SIZE-1];
  346.   int i;
  347.  
  348.   /* Set up the (symmetric tridiagonal) matrix, where the only
  349.    * elements of interest are those on the diagonal. These are
  350.    * stored in array "a[]". Array "d2y[]" initially holds the
  351.    * constants vector, then is overlaid with the calculated
  352.    * variables vector.
  353.    */
  354.  
  355.   m = 6.0/(h*h);
  356.   for(i = 1; i < n-1; i++)
  357.     d2y[i] = m * (y[i+1] - 2.0 * y[i] + y[i-1]);
  358.   a[1] = 4.0 - k_val;
  359.  
  360.   /* Reduce the matrix to upper triangular form */ 
  361.  
  362.   for(i = 2; i < n-2; i++)
  363.   {
  364.     a[i] = 4.0 - 1.0/a[i-1];
  365.     d2y[i] -= d2y[i-1]/a[i-1];
  366.   }
  367.   a[n-2] = 4.0 - k_val - 1.0/a[n-3];
  368.   d2y[n-2] -= d2y[n-3]/a[n-3];
  369.   d2y[n-2] /= a[n-2];
  370.  
  371.   /* Solve using back substitution */
  372.  
  373.   for(i = n-3; i > 0; i--)
  374.     d2y[i] = (d2y[i] - d2y[i+1])/a[i];
  375.  
  376.   /* Solve for end conditions */
  377.  
  378.   d2y[n-1] = d2y[n-2] * k_val;
  379.   d2y[0] = d2y[1] * k_val;
  380. }
  381.  
  382. /* PSPL_COEFF() - Calculate periodic spline coefficients and
  383.  *          return in vector "d2y[]". Matrix to be solved
  384.  *          has the typical form:
  385.  *
  386.  *     +-           -+ +-   -+    +-           -+
  387.  *    | 4  1  0  0  1 | | y1" | = 6 * | y2 - 2*y1 + y0    |
  388.  *    | 1  4  1  0  0 | | y2" |    | y3 - 2*y2 + y1    |
  389.  *    | 0  1  4  1  0 | | y3" |    | y4 - 2*y3 + y2    |
  390.  *    | 0  0  1  4  1 | | y4" |    | y5 - 2*y4 + y3    |
  391.  *    | 1  0  0  1  4 | | y5" |    | y1 - y0 - y5 + y4 |
  392.  *    +-           -+ +-   -+    +-           -+
  393.  *
  394.  *         where m = 6.0/(h*h) and yn" is the second
  395.  *         derivative of the interpolated function at the
  396.  *         "nth" abscissa and y0" = y5".
  397.  */
  398.  
  399. void pspl_coeff(y,d2y,h,n)
  400. float y[];
  401. double d2y[],
  402.        h,
  403.        fabs();
  404. int n;
  405. {
  406.   double c,
  407.      m;
  408.   static double a[MAX_SIZE-1],
  409.         b[MAX_SIZE];
  410.   int i;
  411.  
  412.   /* Check for matching end ordinates */
  413.  
  414.   if(fabs(y[n-1] - y[0]) > 0.0)
  415.     error(NMT_ORD,NULL);
  416.  
  417.   /* Set up the matrix, where array "a[]" holds the diagonal
  418.    * elements, "b[]" holds those elements of column "n-1", and
  419.    * "c" holds the current element of interest of row "n-1".
  420.    * Array "d2y[]" initially holds the constants vector, then is
  421.    * overlaid with the calculated variables vector.
  422.    */
  423.  
  424.   d2y[0] = 0.0;
  425.   m = 6.0/(h*h);
  426.   for(i = 1; i < n-1; i++)
  427.     d2y[i] = m * (y[i+1] - 2.0 * y[i] + y[i-1]);
  428.   d2y[n-1] = m * (y[1] - y[0] - y[n-1] + y[n-2]);
  429.   a[1] = 4.0;
  430.   b[1] = 1.0;
  431.   b[n-2] = 1.0;
  432.   b[n-1] = 4.0;
  433.   c = 1.0;
  434.  
  435.   /* Reduce the matrix to upper triangular form */ 
  436.  
  437.   for(i = 2; i < n-1; i++)
  438.   {
  439.     m = 1/a[i-1];
  440.     a[i] = 4.0 - m;
  441.     b[i] -= b[i-1] * m;
  442.     d2y[i] -= d2y[i-1] * m;
  443.     b[n-1] -= c * m * b[i-1];
  444.     d2y[n-1] -= c * m * d2y[i-1];
  445.     c = -c * m;
  446.   }
  447.   c += 1.0;
  448.   b[n-1] -= c * b[n-2]/a[n-2];
  449.   d2y[n-1] -= c * d2y[n-2]/a[n-2];
  450.  
  451.   /* Solve using back substitution */
  452.  
  453.   d2y[0] = d2y[n-1] /= b[n-1];
  454.   d2y[n-2] = (d2y[n-2] - b[n-2] * d2y[n-1])/a[n-2];
  455.   for(i = n-3; i > 0; i--)
  456.     d2y[i] = (d2y[i] - b[i] * d2y[n-1] - d2y[i+1])/a[i];
  457.  
  458. }
  459.  
  460. /* SPL_INT - Interpolate points using spline function */
  461.  
  462. double spl_int(ix,x,y,d2y,h,i)
  463. float x[],
  464.       y[];
  465. double ix,
  466.        d2y[],
  467.        h;
  468. int i;
  469. {
  470.   double iy,
  471.      t1,
  472.      t2;
  473.  
  474.   t1 = (ix - x[i-1])/h;
  475.   t2 = (x[i] - ix)/h;
  476.   iy = y[i-1] * t2 + y[i] * t1 - SQUARE(h) * (d2y[i-1] * (t2 -
  477.       CUBE(t2)) + d2y[i] * (t1 - CUBE(t1)))/6.0;
  478.   return iy;
  479. }
  480.  
  481. /* IS_FLOAT() - Check that character string is in correct floating
  482.  *        point format. Return TRUE if correct, FALSE
  483.  *        otherwise. The algorithm used is a deterministic
  484.  *        finite state machine. Using the regular
  485.  *        expression terminology of Unix's "lex", the
  486.  *        character string must be of the form:
  487.  *
  488.  *            -?d*.?d*(e|E(\+|-)?d+)?
  489.  *
  490.  *        where d = 0|1|2|3|4|5|6|7|8|9
  491.  */
  492.  
  493. BOOL is_float(str)
  494. char *str;
  495. {
  496.   int c,        /* Next FSM input character */
  497.       state = 0;    /* Current FSM state */
  498.  
  499.   while(c = *str++)
  500.   {
  501.     switch(state)
  502.     {
  503.       case 0:        /* FSM State 0 */
  504.     switch(c)
  505.     {
  506.       case '-':    
  507.         state = 1;
  508.         break;
  509.       case '.':
  510.         state = 3;
  511.         break;
  512.       default:
  513.         if(isdigit(c))
  514.           state = 2;
  515.         else
  516.           return FALSE;
  517.         break;
  518.     }
  519.     break;
  520.       case 1:        /* FSM State 1 */
  521.     switch(c)
  522.     {
  523.       case '.':
  524.         state = 2;
  525.         break;
  526.       default:
  527.         if(isdigit(c))
  528.           state = 2;
  529.         else
  530.           return FALSE;
  531.         break;
  532.     }
  533.     break;
  534.       case 2:        /* FSM State 2 */
  535.     switch(c)
  536.     {
  537.       case '.':
  538.         state = 4;
  539.         break;
  540.       case 'e':
  541.       case 'E':
  542.         state = 5;
  543.         break;
  544.       default:
  545.         if(isdigit(c))
  546.           state = 2;
  547.         else
  548.           return FALSE;
  549.         break;
  550.     }
  551.     break;
  552.       case 3:        /* FSM State 3 */
  553.     if(isdigit(c))
  554.       state = 4;
  555.     else
  556.       return FALSE;
  557.     break;
  558.       case 4:        /* FSM State 4 */
  559.     switch(c)
  560.     {
  561.       case 'e':
  562.       case 'E':
  563.         state = 5;
  564.         break;
  565.       default:
  566.         if(isdigit(c))
  567.           state = 4;
  568.         else
  569.           return FALSE;
  570.         break;
  571.     }
  572.     break;
  573.       case 5:        /* FSM State 5 */
  574.     switch(c)
  575.     {
  576.       case '+':
  577.       case '-':
  578.         state = 6;
  579.         break;
  580.       default:
  581.         if(isdigit(c))
  582.           state = 7;
  583.         else
  584.           return FALSE;
  585.         break;
  586.     }
  587.     break;
  588.       case 6:        /* FSM State 6 */
  589.     if(isdigit(c))
  590.       state = 7;
  591.     else
  592.       return FALSE;
  593.     break;
  594.       case 7:        /* FSM State 7 */
  595.     if(isdigit(c))
  596.       state = 7;
  597.     else
  598.       return FALSE;
  599.     break;
  600.     }
  601.   }
  602.   return TRUE;
  603. }
  604.  
  605. /* ERROR() - Error reporting. Returns an exit status of 2 to the
  606.  *         parent process.
  607.  */
  608.  
  609. void error(n,str)
  610. int n;
  611. char *str;
  612. {
  613.   fprintf(stderr,"\007\n*** ERROR - ");
  614.   switch(n)
  615.   {
  616.     case ILL_ARG:
  617.       fprintf(stderr,"Argument must be greater than zero: %s",
  618.       str);
  619.       break;
  620.     case ILL_CMB:
  621.       fprintf(stderr,"Cannot use -k option with -p option");
  622.       break;
  623.     case ILL_KVL:
  624.       fprintf(stderr,"Value for -k option cannot equal 4.0: %s",
  625.       str);
  626.       break;
  627.     case ILL_OPT:
  628.       fprintf(stderr,"Illegal command line option: %s",str);
  629.       break;
  630.     case INS_INP:
  631.       fprintf(stderr,"Insufficient input data");
  632.       break;
  633.     case MIS_KVL:
  634.       fprintf(stderr,"Missing value for -k option");
  635.       break;
  636.     case MIS_XVL:
  637.       fprintf(stderr,"Missing value for -x option");
  638.       break;
  639.     case MIS_YVL:
  640.       fprintf(stderr,"Missing ordinate value");
  641.       break;
  642.     case NMT_ORD:
  643.       fprintf(stderr,"End ordinates do not match");
  644.       break;
  645.     default:
  646.       break;
  647.   }
  648.   fprintf(stderr," ***\n\nUsage: spline [-aknpx]\n");
  649.   exit(2);
  650. }
  651.  
  652. /*** End of SPLINE.C ***/
  653. 
  654.   fprintf(stderr,"\007\n*** ERROR - ");
  655.   switch(n)
  656.   {
  657.     case ILL_ARG:
  658.       fprintf(stderr,